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We study the Richardson equations close to the critical values of the paring strength g c where the 
occurrence of divergencies preclude numerical solutions. We derive a set of equations for determining 
the critical g values and the non-collapsing pair energies. Studying the behavior of the solutions 
close to the critical points, we develop a procedure to solve numerically the Richardson equations 
for arbitrary coupling strength. 
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O ' I. INTRODUCTION 

o : 

The exact solution of the BCS or pairing Hamiltonian was presented by Richardson in series of papers beginning in 
Qh' 1963 [lj,|2f : just a few years after the seminal paper of Bardeen, Cooper and Schrieffer [3J. Despite of the new avenues 
of research that it could have opened at this early stage in the development of the theory of superconductivity, the 
work passed almost unnoticed due, perhaps, to the technical difficulties in solving numerically the set of non linear 
equations for the spectral parameters. In recent year, the Richardson exact solution was rediscovered and applied 
successfully to ultrasmall metallic grains where it was shown to be essential for the description the soft crossover 
between superconductivity and the paring fluctuation regime as a function of the grain size Since then the 
Richardson model has been extended to a wide class exactly solvable models called Richardson-Gaudin (RG) models 
^v^j ' with potential applications to several quantum many body fermion and boson systems 0, 0, Q> Q- However, the 
numerical treatment of the exact solution for moderate to large size systems is still a cumbersome task in spite of the 
Q\ • recent efforts to overcome this problem. 

For the sake of clarity, let us begin by introducing the Richardson equations for an M fermion pair system: 

\Q' 

•S l-4 5 £-^L + 4 ff £ =0, (1) 

J* (*¥a) 

where g is the pairing strength, rjj are a set of L parameters usually related to the single particle energies, dj are 
the effective degeneracies defined as dj = fj/2 — with Vj de number of unpaired fermions in level j (Seniority 
quantum number) and flj the total degeneracy of the level j. The e a 's are the M unknowns parameters called pair 
energies. Given the L parameters rjj, the L effective degeneracies dj and a pairing strength g, the pair energies e a are 
obtained by solving the set of M nonlinear equations (QJ. However, it is in general not easy to find a good initial guess 
that would lead directly to the appropriate solution. Instead, it is customary to begin with a given configuration in 
the weak coupling limit (<? — > 0) where the equations Q can only be satisfied for e Q — * 2r)j. The exact solution is then 
evolved step by step for increasing values of g up to the desired value [j| . In each step one uses the solution obtained 
in the previous step as the input data for the numerical solver. Even if the procedure is successful, leading to the 
correct solution, it involves a heavy numerical work. In most cases, this procedures is stopped due to the existence 
of singularities for some critical values g c of the pairing strength. At g — g c a subset of pair energies e a turn out to 
be equal to 2rjk for some k, giving rise to divergencies in some terms of eq. Q |lOj |- Remarkably, these divergencies 
cancel out and the corresponding solution (the set of e a ) is a continuous function of g in a neighborhood of g c . In 
fact, the solutions are always continuous for every value of g. 

The existence of such critical points has significant consequences in the numerical solution of the equations which 
becomes unstable near g c . Very slow convergency, no convergency at all or jumping to another solution are the typical 
problems found close to g c . Moreover, there may be more than just one critical value in the real interval [0, g], and we 
do not know a priori where the g c 's are located. Those values are known only after one has already solved the equations 
for a large number of points surrounding the <ji c 's and making some kind of interpolation afterward. Therefore, the 
localization of the critical values of g relies on heuristic procedures and is not based on any mathematical properties 
of the equations. 

An alternative procedure to cross the critical region, based on a non linear transformation of the collapsing pair 
energies was recently proposed 11]. However, this procedure is unable to predict the critical values of g. 

In this paper we will analyze the properties of Richardson equations in the vicinity of critical g values. We will 
derive conditions that allow to a priori determine all the critical values g c associated to any "single particle level" rjj , 
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and provide the exact solution at these points. In addition, we will describe the asymptotic behavior of the solution 
in the limit g —* g c , and we will present an algorithm to solve the Richardson equations for values of g near any g c . 



II. TRANSFORMING THE RICHARDSON EQUATIONS. THE CLUSTER EQUATIONS 

Our approach to deal with the singularities in eq. JJ| consists in transfor ming the system through a change of 
variables first suggested by Richardson [iJJ] and later used by Rombouts et al. [TI| to develop a numerical algorithm 
for solving the equations near the critical points. We have already mentioned that for critical values of g some subset 
of the pair energies becomes equal to one of the values of 2r)k- It can be shown that the number of such pair 
energies is Mk = 1 — 2dk- We can then characterize a critical point g c by the condition lim s ^ Sc e a — 2r\k Va £ Ck, 
where Ck stands for the subset of indices of the Mk pair energies that satisfy that limit, i.e. the e Q 's that cluster 
around the real point 2r\k in the complex plane for g near g c (see |l2j for a graphical representation of these clusters). 
Therefore, we will deal with two subsets of variables, the Mj~ e Q 's with a £ Ck that give rise to the singularities, 
and the remaining (M — Mk) variables with a ^ Ck which behave smoothly close to g c . Consequently, we will treat 
separately the Mk Richardson equations with a G Ck (the cluster equations): 

l-Ag-Jt—Ag ]T £ -^+4.g£ -J— =0. (2) 

2 Vk -e a ^ 2 Vj - e a ^ e a - ep ^ e a - ep 

The second and fourth terms of eq. (j2J diverge for g — > g c since (2r]k — e a ) and (e a — ep) go to zero. Moreover these 
quantities must approach zero at the same rate in order to cancel out. To avoid the singularities we will multiply the 
cluster equations by (2rjk — e a ) p (for some p > 0). With this in mind we will introduce a change of variables for the 
pair energies in the cluster (see ref. 0>E1) 

S p = ( 2 »ft " e -) P for P = 1, 2, . . . , M k . (3) 

aec k 

This is an invertible transformation and, in principle, we can always recover the e a 's for any arbitrary set of S p . 
Keeping in mind that we have just Mk independent variables S p we will extend the definition allowing p to be any 
positive integer or zero (note that So = Mk). The variables S p behave smoothly in the vicinity of g c , they are real 
and, for p > 0, they satisfy lim g ^ gc S p = 0. 

In order to transform the cluster equations @ we first multiply them by (2rjk — e a ) p , for a generic integer power 
p, and then sum the resulting equations over a S Ck- 

S p+ 2 9 pS p ^-2i£s p ^S i -A g £ £ rfj( f fc - e " )P +4 g £ £ (2 ^ fc = Ca)P = 0. (4) 



Note that the last two terms in the left hand side of the equation cannot be easily rewritten in terms of S p if we 
restrict to values of p < Mk- To overcome this limitation, we will make a series expansion of those terms, allowing p 
to be any integer. We are interested in the limiting behavior of the solution close to g c so that our reasoning will be 
valid in an interval g £ (g c — 8, g c + 8) for some small enough radius 8 > 0. Taking into account that (2r\k — e a ) is 
infinitesimal at g c and that ep in equation (0} does not belong to the cluster, we have: 

1 1 ^ / 2i] k - e a 



2r]j-e a 2i] 3 ■. - 2r lk ^ \2rj k - 2-q 

1 = 1 y> / 2r] k - e a \ " 
e a - ep 2r) k - ep \ 2r\ k - ep J 

These are absolutely (and rapidly) convergent geometric series, a property that will remain valid after a finite 
summation over a. Introducing eq. JSJ in eq. Q and making the summation over a we get for p = 1: 

oo 

Sx - 4 3 d k M k - 2g(M 2 k - M k ) + 4g ^ S 1+n P n = 0, (6) 
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where the {P n } is a set of coefficients that will be defined below. Taking in © the limit g — > g c the second and third 
terms must cancel out and one gets the cluster condition Mk = f — 2d k . Similarly, for p > 1 we obtain: 

p— 2 oo 

S p - 2g(M k + 1 - p)S p - 1 - 2. 9 ^ S p _i_iSi + 4 5 S P+ „P„ = for all p > 1. (7) 

i=l ra=0 

The coefficients P n appearing in © and J7J are defined as: 

The infinite recurrence relation (J7J| is the fundamental equation of our formalism. In order to achieve our main goals 
we need to establish some important properties (lemma ^ and lemma 0) of the variables S p and the series appearing 
in eqs. © and 10. It is easy to show that these series are bounded by their first element S p through the condition: 



^ S p +nPn 



n=0 



< \SJ D ■ 
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(9) 



where D is a positive bounded number involving rjj and ep, and x is an upper bound for the absolute values of the 
geometric series ratios in eq.© satisfying lim g ^ gc x — 0. On the other hand, we have that \S p \ < J2aec k |2'7fe — e a\ p 
with (2i]k — e a ) — > as g —> g c . Therefore, we will assume in what follows that there is an integer number n$ big 
enough such that we may disregard every term S p (and related series) for p > ng. 

The next aspect we need to analyze is the order of the infinitesimal S p in relation to the order of Si. We will 
demonstrate that Si, S2, S3, . . . , SM k are of the same order, implying that lim 9 ^ 9c S p /Si = Xp f° r P !- Mk- On 
the other hand, S p is of higher order than Si if p > Mk, that is lim ff ^ ffc S p /Si — 0. In lemmas ^ and |3 we will 
demonstrate both assertions. 

Lemma 1. Assuming that (1 + 4gPo) ^ 0, with g e (g c — 5, g c + S) and Pq defined in 0), the set {S2, S3, . . . , S M } 7 
for some positive integer [i with 2 < /z < Mk, are infinitesimal of the same order as Si at g = g c . In the general case 
fi = M k . 



Proof. Rearranging equation Q we get 

1 



S n 



(1 + 4.9P0) 



p-2 



2g{M k + 1 - p)S p -i + 2gJ2S P -i-iSi - 4 5 S p+n P n 



(10) 



4—1 n— 1 

To prove our statement we will apply this equation recursively starting with p = 2. In this case we get: 

_ 2g(M k -l) 4 ff 



-Si 



(1+4 5 P ) (1 + 4 9 P )^ 



5> 



2+nPn, 



so that S2 is of the same order as Si. In the next step we set p = 3 in eq. H10|l . getting S3 in terms of S2 and S p 's 
with p > 3. Replacing S2 with the value obtained previously and rearranging the resulting expression we get for S3: 



S 3 



1 



(l + 4gPo) 2 + 8g 2 Pi(Mk-2] 



[(2g) 2 (M k - 2){M k - l)Si + 2g{\ + 4gP )S 2 



-4g VJ(2 5 (M fe - 2)P n+ i + (1 + 4.gP )P„)S 



3+n 



from where we see that S3 is of the same order as Si. We can continue this process replacing in each step S p _i in eq. 
I|10(l for the expression computed in the previous step and solving again for S p , getting for every p > 2 and p < Mk'. 



S p = B p Si+J2 D n )S P+n+ J2 X S S * S . 



(11) 



n=l 



i+j<p 
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where the coefficients B p , D^j/ and Xf) can be calculated recursively. For intance, for the linear terms we have the 
recurrence relation 



(1 + AgP ) - 2g(M k - p)D™ 



D 



(p+i) 



2fj 



(l + AgP Q )-2g(M k -p)D [ { 



(p) 



((M k ~p)D { n P l 1 -2P n ), (12) 



starting with the known coefficiens for p = 2. 

In this way we have shown that the local expansion of S p has linear terms in Si and in S p + n plus quadratic terms, 
thus we conclude that S p is an infinitesimal of the same order as Si for 1 < p < M k . 

Let us now consider what happens when p takes the value M k + 1. In such a case the factor (M k + 1 — p) in eq. 
becomes zero and the chain of replacement is broken. The expansion of SV fe +i has only linear terms of higher order 
and quadratic terms involving Si . Therefore for p > M k there are no linear terms in Si . We could have reached the 
same conclusion analyzing the recurrence relation l |12|l . 

Finally and for the sake of completeness, let us mention that after replacing S p ~i in equation (|10|l the terms involving 
S p on both sides of the equation might cancel out. It would be equivalent to the cancellation of the denominators in 
the recurrence relation (|12H . In such a case our process of replacement is interrupted and ends for p = p < M k . □ 

So far, we have shown that the first S p s are of the same order as Si. At the same time, is was suggested that S p 
is of higher order if p > M k . In the next lemma we will proof this statement. 

Lemma 2. The set of variables {S p } for p > M k are infinitesimal of higher order than Si at g = g c . 

Proof. In this case we proceed in the reverse way as we have done in lemma ^ We will begin with an S p for a large 
value of p and carry out the replacement in backward direction. Once again we start rewriting equation Q to obtain: 



Sn 



P (M k -p) 



4 p \ 00 p— i 

$p+i + 2 Sp+i+ n P n — S p 



'Mi 



n—1 i—1 



(13) 



We will now take, for g € (g c — S,g c + 8), an integer number n$ large enough such that we may disregard every 
term S p for p > n$. We start making p — n$ in eq. (|13fl . After removing the negligible terms it leads us to: 

^ ns-l 

Sns = (M k -n 5 ) ^ Sn *- iSi - 

In the second step we make p = rig — 1 in eq. (| 1 31) and replace the value of S ns computed before. Repeating the 
same operation for p = ns — 2, rig — 3, . . . , M k + 1 and replacing, in each case, the linear term for its value computed 
in the previous step, eventually we get the following quadratic expression for every p > M k : 



i+3>P 

(p) 

In this equation the coefficients Y^j can be calculated recursively and the indices i, j take values starting from 
1 but keeping the condition i + j > p. As only quadratic terms are retained we conclude that S p for p > M k are 
infinitesimal of higher order than Si , S2, ■ ■ ■ , SM k ■ Note once again that the chain of replacement cannot be continued 
for p < M k . □ 

So far we have shown some essential properties of the variables S p and the related series that are required to develop 
our formalism. In the next section we will derive the equations that allow to calculate the values of g c associated with 
any single particle level rj k . 



III. CRITICAL VALUES OF g AND THE SOLUTION OF RICHARDSON EQUATIONS 

In this section we will derive a set of equations suitable to compute all the values of g c for any single particle level 
r\ k . Furthermore, we will obtain, at the same time, the solution of Richardson equation at g c , namely, the values of 
the unknowns ep for fj ^ C k , assuming that the variables in the cluster reach their limiting value e a = 2 r\ k (e a G C k ). 
We summarize our results in the following theorem: 



5 



Theorem 1. All critical g values of the Richardson Equations 0) associated with the single particle level and 
the corresponding values of the non collapsing pair energies ep ((3 C k ) are the solutions of the following system of 
equations: 



(1 + 4 5 P ) igPi 
-2g{M k -\) (1 + 4 5 P ) 
-2g(M k - 



AgP 2 
4gPi 

(1 + 4.gPo) 



• ' 4 5 P Mfe -i 

• • 4 3 P Mfc -2 

•■ 4<?P Mfc - 3 

2g (1 + 4 5 P ) 



= 



(15a) 



d 1 1 

l"4gV^ +4gM k — + 4g V = for all a $ C k (15b) 

^-J 2r)j -e a e a - 2r] k e a - ep 

(Pita) 

Proof. Let us consider the transformed Richardson equations according to the change of variables denned in eq. 10 . 
This system is formed by the (M — M k ) Richardson equations P for a ^ C k together with eq. © and the subset 
of equations for p = 2, 3, . . . , M k . It is a set of M independent equations equivalent to the original ones. In fact, 
we can consider the unknowns Si, S2, ■ ■ ■ , SM k as independent variables together with {ep} for /? ^ C' k . The variables 
{e a } for a € C k , that still appear in the equations, and the S p for p > M k are formally functions of Si, S2, ■ ■ ■ , SM k 
computable through the inverse of transformation Q . At this point we are ready to analyze the M k cluster equations 
(0 and [7] for p < M k ) in the limit g — > g c , keeping the lower order terms. According to lemmas and J5J we can 
discard terms with p > M k and terms involving the products Si Sj , since they are of higher order than Si , S2 , ■ ■ • , SM k , 
thus we get the system of equations: 

(1 + AgP )Si + AgPiS 2 + AgP 2 S 3 + ■■■+ igP Mk -iS Mk = 
-2g(M k - l)Si + (1 + 4gP )S 2 + 4gPiS 3 + ■■■+ AgP Mk _ 2 S Mk = 

-2g(M k - 2)S 2 + (1 + 4 5 P )5 3 + • • • + AgP Mk ^S Mk = (16) 

-2gS Mk -i + (l + 4,gPo)S Mk = 

This is a linear homogeneous system with the trivial solution Si = 0, S 2 = 0, . . . , SM k = 0, and indeed the variables 
S p vanish for g = g c . Since we are looking for nontrivial solutions that are continuous functions of g valid also in 
the vicinity of g c where S p ^ 0, the condition for such solutions to exist is the cancellation of the determinant of the 
linear system (shown in eq. 115b). As the coefficients of the system are functions of the variables e Q not belonging to 
the cluster we have to resort to the (M — M k ) remaining Richardson equations with a ^ C k . After realizing the 
limit g — > g c , we obtain the system of equations 115fl introduced by the theorem. □ 

As we have already stated, the numerical solutions of the non-linear system of equations (jl5|l provide all the critical 
values g c for the single particle level r\ k that has been previously selected. In addition, for each solution, we get the 
complete set of pairs energies e a denning the wave function of the quantum many body system. 

With respect to the possible complex g c solutions, they may be still interpreted as critical values for not hermitian 
integrable hamiltonians, but we will not analyze these cases in the present paper. 

Once we have the numerical solution of Eqs. (|15J) we can use the values of g c and e Q 's (a ^ C k ) to compute the 
coefficients of the homogenous linear system l|16f) . Solving this linear system yields (after normalizing to Si = 1) the 
limit ratio \v = um 9^9c sf- As S p = at g = g c for all p > 1, we can obtain from eqs. (|16[) the derivatives: 

(17) 

*1 n = n„ 



IV. SOLVING RICHARDSON EQUATIONS NEAR g c 



So far we have shown how to compute the solution at g = g c . We will now derive a method to approach the solution 
at a value go close to the critical point g c . The main issue in the numerical solution of the Richardson equations 
is to determine a good initial guess for the pair energies, specially in the vicinity of g c where the equations become 
unstable. With that initial guess the solution at g — go is obtained with standard numerical technics. Once we have 
the solution for this specific value go, we can reach a more distant value of g, by increasing (or decreasing) g, step by 
step, using the solution of the previous step as the starting guess for the next one. 
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An appropriate set of starting values at go = g c + Sg can be obtained by means of a linear approximation for the 
parameters S p and e a : 

Sl ^o)-(^^j Sg, S p (go) « *p ^ (p > 1), 



e Q (ffo) ~ e a (g c ) + 



de a 
dg 



Sg (a <£ C k ), 



(18) 



9=9c 



this approximate values will be the initial guess to solve the Richardson equations. 

In order to determine the derivative dSi/dg at g = g c we will consider up to second order terms in equations © 
and J7J and substitute the variables S p in terms of the quadratic expression S p ~ XpSi + a-pSf, where \p = f° r 
p > Affe, xi = 1 and ai = 0. From equation (|14|l we can see that S p has quadratic terms in Si only for p < 2Mk, 
while for greater values of p the order in Si is higher due to the condition i + j > p. Therefore, we will retain just 
the first 2M k variables S p and deal with a system of 2M k equations. Computing the derivative with respect to g of 
the resulting system, and after making some simple algebraic manipulations, we get for p = 1: 



2M k 



71=1 



and for 1 < p < 2M k 



p-2 



a p - 2g(M k + 1 - p)a p -i - 2g ^(Xp-i-i Xi) 



2M k -p 



4 5 l(Xn+p + a-n+p Si)P' n + a n+p P n S-(} 



X P + a p Si 



71=0 



where the primes stand for the derivatives with respect to g and, according to definition JSJ, we have: 

, = y> (n+ 1) / dep' 

Taking the limit g — > g c and reordering the equations in a convenient way we obtain 

/ , M k \ 2M k 



for p = 1, and 



^ - 2 5c S/ Efe-i-i Xi) + 4. 9c ^ Xn+p-P,; ) - 2g c (M k + l-p) a p ^S( 



M k -p 



(19) 



(20) 



(21) 



(22) 



2M k 



(1 + 4g c P ) a p S; + ]T {Ag c P^ p ){ aj S{) = 0, 
j=p+i 



(23) 



for 1 < p < 2M k . This is a non-homogeneous system of 2M k equations with a set of 2M k + (M — M k ) unknowns: S*/, 
tt2, 03, • ■ ■ , a,2M k joined with the derivatives {dep/dg), for [3 <^C k . 

The system of equations (12. '{t is non-linear because of the products ajS{. However, we can obtain a linear system 
for the derivatives with a unique solution by defining a matrix B as: 



/ P-2 M k -p \ 

B p,l = ( ~ ~ 29cSi^2(Xp-i-l Xi) + E Xn+pPn , 

\ 9c i=l n=0 / 



B PiP _i =-2g c (M fc + l-p), 
B p j = Ag c Pj-p {j > p), 



B P , P = (l + 4 9c P ) (p>l), 
S P)j = (Kj<p-1); 
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and the non-null vector v : 

dSi dSi dS\ 

1, a 2 — , a 3 — , . . . , a 2 M k -j- 

dg dg dg 

Using these definitions, equations (|23[) can be rewritten as: 

B v = 0, 

such that the following condition must be satisfied: 

dct(B) = 0. (24) 

Since the derivatives appear just in the first column of the matrix B, (|24|l is a linear equation in ^ and in the 

(M — Alk) derivatives ^j- (for (3 ^ C k ) that does not include the unknown coefficients a p . A set of (M — M/.) equations 
is required to complete the system. In order to get the necessary equations we will compute the derivatives of the 
Richardson equations out of the cluster: 

l-4^-^+4 9 ^— +4^— = {aiC k ). (25) 

U 2 ^~ ea Hc k ea ~ ep pTi k ea ~ ep 

Expanding the last term in the variables S p we obtain 

L , oo 

Computing the derivative with respect to g and taking the limit g — > g c we finally get: 



1 L 

17. t—^ 




1 f dep de a 



"^TT^-O { *iC k ). (27) 



(2% - e Q ) 2 V dg J g=gc ^ (2^ - e a )^ J dg 

Equations (|24|l and H27J) constitute a non-homogenous linear system with a unique solution for the desire derivatives. 

V. A NUMERICAL EXAMPLE: FERMIONS IN A 2-DIMENSIONAL LATTICE 

In this section we will apply our formalism to a pairing model of fermions in a 2-D square lattice of N x N sites with 
periodic boundary conditions. This example was previously treated in Ref. for the ground state with a repulsive 
pairing interaction (g > 0). We will now analyze the model using the methodology developed here, focussing on the 
critical g values and in the behavior of the solutions in their vicinities. 

We will consider a 6 x 6 lattice with 36 fermions (half filling, M — 18). The pairing Hamiltonian in momentum 
space is 

h p = E £ka k ak + f E 44 a k' ak ' 

k k 

k' 

where k = [k x , k y ] and £k = — 2[cos(fe c ) + cos(k y )) = rfc, and k is the time reversal of k. The single fermion energies 
Ej and the corresponding degeneracies flj for the 6x6 lattice are displayed in the Table [I] In the limit g — the 
ground state is obtained by distributing the M = 18 pairs in the lowest possible states. In this numerical example 
we will work within the fully paired subspace containing the ground state ( vj =0, dj — —Clj/4 ). Consequently, the 
excited states to which we will refer to in this example are within the seniority subspace. 
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FIG. 1: Real part of the pair energies e a for attractive (g < 0) and repulsive (g > 0) pairing in a 6 x 6 lattice. 
TABLE I: Single fermion energies and degeneracies for the 6x6 lattice. 
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We will begin solving the system of equations ljl5JI for a specific single particle level, for instance, j = 4 and £4 = — 1 
(see table HJ. The degeneracy of this level is 8 and the number of pair energies in the cluster that collapses at 2e4 — — 2 
is given by the cluster condition M4 = 1 — 2d,4 = 5, one more than the value allowed by the Pauli exclusion principle. 
The solution of the system (|15|l provides all the information concerning the many body state at each specific critical 
value of g. In table[n]we show, for j — 4, the first negative g c 's together with the corresponding energy eigenvalues. 

A different kind of analysis is presented in Table ITTT1 where we show the first critical g values for the ground state. 
In this case the clusters collapse at different single particle levels for each g c as indicated. All the cluster are formed 
by 5 pair energies except for j = 1 where Mi = 3. In order to have a global vision of the results we show in Fig. Q]the 
real part of the pair energies for the ground state solution of the Richardson equations Q for positive and negative 
g values. The full lines are the pair energies e a and the horizontal dotted lines correspond to 2s j. The critical points 
of Table ITTT1 can be seen in the figure at the crossing point of each cluster with twice the corresponding single particle 
energy. 

The solutions near the critical points g c were obtained following the approach described in Section llVl With the 
solution for the first g c next to g = (negative or positive depending on the case), i.e, the ep for (3 ^ Ck and value 
of g c itself, we solve the linear system of equations l|24[) and l|27[l to obtain the derivatives (dSi/dg) and (dep/dg) at 
the critical point. In addition we solve the linear system to get the coefficients \v With all this information 
we determine a good initial guess at some go near g c by mean of the linear approximation l|18(l . Next we solve the 
Richardson equations at go, a process that rapidly converge to an accurate numerical solution. With the solution at 
go as a starting point we move to the next value of g = go + Ag. Proceeding in this way and updating the starting 
values at each step we compute all the desired points in the curves. We repeat the process for each critical g. The 
pair energies obtained from the solutions for g « smoothly connect with those obtained near the first g c . The same 
smooth behavior is observed between two consecutive critical values. 

The good quality of the initial guess near g c is due to the smooth behavior of the variables S p and ep (/? ^ C'k) that 
allows us to use the linear approximation l|18|) . For the e^'s this fact can be appreciated in figure ^ In figure [5] we 
show a graphical representation of the first S p 's in a range of g around g c — 0.170878 corresponding to the collapse 
at the single particle level j = 2, 2e2 = —6 of a cluster of Mi = 5 pair energies. We can we see in the figure the 
linearity of the S^'s for 1 < p < 5 confirming that S p is an infinitesimal of the same order as Si for 2 < p < M2. In 
addition, we have plotted the results for Sg, which clearly shows that it is an infinitesimal of a higher order. 



0.15 0.16 0.17 0.18 0.19 0.20 

g 



FIG. 2: Variables S p close to the collapse at g c = 0.170878. The first 5 S p 's are displayed in solid lines while Se is drawn in 
dashed line. 



TABLE II: First critical g values for the cluster collapsing at 2e4 = —2 





Energy 


-0.0384565 


-47.6184 


-0.0391412 


-49.5405 


-0.0394719 


-53.3549 


-0.0404240 


-55.5262 


-0.0412922 


-57.4106 


-0.0413245 


-62.5795 



VI. CONCLUSIONS 

In this work we have studied the solution of the Richardson equations close to the critical values of the coupling 
constant g. We have derived a set of well behaved equations to determine the actual critical values g c associated 
with any single particle energy , and the complete set of pair energies defining the exact eigenstate of the system. 
In addition, we studied the behavior of the solutions close to the critical points, obtaining a linear approximation 
for the pair energies that serves as a good initial guess in the critical region. With this formalism, one can solve 
numerically the Richardson equations around the critical points, overcoming the numerical instabilities that usually 
arise. The knowledge of the critical values of the coupling constant greatly simplifies the numerical treatment of the 
equations. To illustrate the formalism we have analyzed the pairing Hamiltonian in a 6 x 6 square lattice at half filling 
for repulsive and attractive pairing strength. Making use of our approach we have located the g c 's for the ground 
state and several excited states in the Seniority subspace. The numerical solution between consecutive values of 
f/ c 's can be easily carried out by solving the original Richardson equations (Q. With this new approach we hope to 
be able to solve exactly large systems with arbitrary single particle energies and degeneracies. 



TABLE III: Critical g values for the ground state 





2e 4 = -2 


2e 3 = -4 


2e 2 = -6 


2ei = -8 


g c (negative) 
g c (positive) 


-0.0413245 
0.598232 


-0.0635021 
0.240579 


-0.0877434 
0.170878 


-0.131927 
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